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Abstract. "Coil design" is an inverse problem in which arrangements of wire are 
designed to generate a prescribed magnetic field when energised with electric current. 
The design of gradient and shim coils for magnetic resonance imaging (MRI) are 
important examples of coil design. The magnetic fields that these coils generate are 
usually required to be both strong and accurate. Other electromagnetic properties of 
the coils, such as inductance, may be considered in the design process which becomes 
an optimisation problem. The maximum current density is additionally optimised in 
this work and the resultant coils are investigated for performance and practicality. 
Coils with minimax current density were found to exhibit maximally spread wires 
and may help disperse localized regions of Joule heating. They also produce the 
highest possible magnetic field strength per unit current for any given surface and 
wire size. Three different flavours of boundary element method that employ different 
basis functions (triangular elements with uniform current, cylindrical elements with 
sinusoidal current and conic section elements with sinusoidal-uniform current) were 
used with this approach to illustrate its generality. 
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1. Introduction 

In magnetic resonance imaging (MRI) an exquisitely uniform and very intense magnetic 
field is used to polarize the spin population of a sample so as to maximise the strength 
of the nuclear magnetic resonance (NMR) signal. A process known as "shimming" is 
performed at the start of every scan to ensure that this field is as uniform as possible. 
Shimming involves adjusting the electric current in a set of "shim coils" that each 
generate a magnetic field of spherical harmonic intensity in the region of interest (ROI). 
MR images are formed by superimposing magnetic field gradients which causes the 
frequency of the NMR signal, the Larmor frequency, to vary linearly across the sample. 
Fourier techniques reconstruct the image from these frequency encoded NMR signals. 
The linearly varying magnetic fields are generated by "gradient coils" . This paper deals 
with the design of both gradient and shim coils that dictate the speed, resolution and 
accuracy of MRI [1] . 

Gradient and shim coil design is an inverse problem in which arrangements of 
wire are required to generate a specified magnetic field when energized. Additional 
considerations are required such as minimal stored energy, so that they may be switched 
rapidly, or minimal resistive power dissipation, so that their temperature does not 
increase excessively. As MRI machines get shorter to improve patient comfort [2] so 
too must the gradient coils [3]. Reducing the length of the gradient coils pushes the 
wires closer together to maintain magnetic field accuracy. However, gradient and shim 
coils are constructed from finite sized wire and so there is a minimum wire separation 
that can be built. In this work the maximum current density was minimized in the coil 
design process which maximally increases the minimum wire spacing of a coil for fixed 
coil surface geometry. For a given engineering limit for the minimum spacing between 
wires this technique can be used to increase the efficiency of the coil (the amount of 
field per Ampere). It can also be used to reduce the local power dissipation and disperse 
the hot spots of a coil. The present study demonstrates the design of coils with some 
increase in inductance or resistance in order to spread wires. Such designs should be 
judged by appropriate metrics that better encapsulate the coil design problem than 
those designed to reflect purely the stored energy or power dissipation of the design. 

Coil design is also known as magnetic field synthesis and is described by a Fredholm 
equation of the first kind, which is known to be ill-posed [4j. Early approaches to coil 
design in MRI cancelled undesired spherical harmonic components of the magnetic field 
by symmetry and appropriate positioning of loops and arcs of wire (e.g. [5]) or by 
parameterized surface current densities [6]. The "target field" method [7] circumvented 
the problem of ill-posedness by employing a Fourier-Bessel expansion of the | r ^ r , Green's 
function (see §3.11 in Ref. [8]) and defining a continuous target field function for 
— oo < z < oo. Since Fourier transforms have unique inverses, it is possible to 
analytically determine the current density on an infinitely-long cylinder for a limited 
set of target field functions. 

If the target field is defined at a finite set of points, there exists an infinite number of 
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current densities that can produce such a field. With the introduction of the minimum 
inductance as a constraint [9] the problem becomes regularized and a unique solution 
exists. This is essentially Tikhonov regularization [10J of the ill-conditioned system 
of linear equations [HI [121 EE]- This minimum inductance solution was found to be 
somewhat impractical, so the field was allowed to deviate from its prescribed target 
values in order to permit a smoother wire pattern. In a related method, a smooth current 
density was always obtained if it was defined as a weighted sum of a finite number of 
truncated sinusoidal functions [El [151 EJ [El EE] . Parameterizing the current density 
in this manner allowed more practical, finite-length cylindrical coils with limited spatial 
frequency to be designed in a method sometimes referred to as the Turner-Carlson 
method. Minimization of the inductance (the current-normalized stored energy) can be 
easily substituted by the resistance (the current-normalized power dissipation) while still 
resulting in a unique solution [16] . For coil designs with asymmetric target field location 
it is necessary to enforce zero net torque in the presence of an intense background 
magnetic field [17] . 

A coil may be defined by its surface current density. The magnitude of that current 
density defines the wire spacing and the amount of heat generated at positions on the 
surface. To spread the closest wires or to reduce the local heating, the maximum value 
of the current density magnitude was incorporated into the coil design problem and its 
maximum value was minimised. This term is not linear nor quadratic, but only convex 
with respect to the current density. In contrast to more traditional approaches where 
only linear systems are solved, we used techniques of convex programming to handle 
the non-linearities and singularities that arise from the "max" term. We developed an 
original algorithm to solve this optimisation problem that can be seen as a continuation 
of two works by Y. Nesterov [El [19]. It can be shown to converge to the global minimizer 
of the cost function but details of this algorithm will be presented elsewhere. 

The concept of minimum maximum current density (minimax|j|) coil design is 
general and is not limited to any particular coil design method. In this work, three 
different boundary element methods (BEM) were used to investigate the behaviour 
of minimax|j| coils. These are the Turner-Carlson [HI [15] . triangular [20] and 
axisymmetric [2T1 [22] BEMs. The target field was specified at a finite set of discrete 
points in a region of interest (ROI). The field synthesis problem is defined as minimizing 
the sum-of-squares field error and is an ill-posed problem. Therefore, a regularizing term 
must be included to obtain a unique solution. 

In a previous attempt to reduce the maximum current density the regularization 
term was adaptively modified [23]. This method showed a considerable reduction in 
maximum current density, but it was not known how optimal the solutions were. Other 
approaches in which regions of the coil were designed manually have been used to control 
the maximum current density: for example, by predefining the return conductors [241 125] 
or by manually introducing a large number of constraints (page 146 of Ref. [21]). The 
method presented here truly minimizes the maximum current density for coils designed 
on surfaces of arbitrary shape that generate any physically realizable magnetic field. 
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2. Methods 

2.1. Physical Model 

In magnetostatics, Ampere's Law, V x B(r) = /i J(r), relates the magnetic field, B(r), 
and the free current density, J(r). Current density must be conserved, so V • J(r) = 0. 
Employing the magnetic vector potential, A(r), where B(r) = VxA(r), Ampere's Law 
becomes a Poisson equation, V 2 A(r) = yU J(r), which has the solution 

A(r) = f [*¥UV (1) 
47r J |r - r'| 

in the Coulomb gauge (V • A(r) = 0). /i is the permeability of free-space and has the 
value 47r x 10~ 7 Hm -1 . With some algebra, this leads to the familiar volumetric integral 
form of the Biot-Savart law [8], 

B(r) = tJ J(r/) x v^f dv '- (2) 

For each directional component of B(r), is a Fredholm equation of the first 
kind which is known to be ill-posed. MRI conventionally requires a strong magnetic 
field for polarization of the nuclear spin states within the sample to be imaged pQ. 
We consider a system immersed in a background magnetic field, Bo(r), that is highly 
uniform, unidirectional and very strong, i.e. Bo(r) = Bq z z, where z is the unit vector 
parallel to the z-axis. The magnitude of the combined magnetic field, \Bq z z + B(r)| ~ 
(B 0z + B z (r)), dictates the local Larmor frequency of the NMR signal and subsequent 
spatial localization. Therefore we only need to design coils to generate a specific B z (r), 
justifiably neglecting the two other components, B x (r) and B y (r). 

In the context of this paper, "coil design" is the inversion of to design an 
arrangement of wires that, when energized, form a current density which generates a 
prescribed magnetic field. The region of space in which the field is prescribed, the ROI, 
is separate from the region in which the current density exists. 

Coil design is rarely as simple as inverting (EJ) but requires the consideration of 
other electromagnetic properties. The stored energy, W, associated with J(r) is [8J 

w = £_[ [ ^ ■ Y W ', (3) 



in 



nJn r 



where Q c is the region of the coil in which J(r) is confined to flow. The resistive power 
dissipation, P, is 

P = pcu [ |J(r)|W, (4) 

Jn c 

where pc u is the resistivity of the conducting medium which, in this case, is assumed to 
be copper, p Cu = 1.68 x 10~ 10 f2 m. 

The coil may be in close proximity to other conducting surfaces defined by the region 
Q e . Changing J(r) in time causes B(r) to also change, Faraday's Law, V x E(r) = 
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and J(r) = crE(r) show that currents may be induced in other conducting surfaces. 
These "eddy currents" can cause deleterious effects on MRI. So, a coil designer must 
consider the effects that the induced eddy currents have on the field in the ROI [26] . For 
low frequencies (<10 kHz) the quasistatic approximation may be used and following the 
approach of Peeren [27] , a Heaviside function response in the coil current was assumed. 
This leads to a linear relationship between coil currents and eddy currents. 

Lorentz forces act on the coil when immersed in a background magnetic field, B . 
The net Lorentz force is zero for a divergence- free current density in a uniform Bo, but 
there may exist a consequential net torque, r, 

r = / r x [J(r) x B ] alV. (5) 
Jn c 

Current density was confined to flow on thin surfaces so that a scalar stream- 
function, VK r )> can be used to define the vector current density 

J(r) = V x [^(r)n(r)] , (6) 

where n(r) in the unit vector normal to the surface at r. 

2.2. Discrete Formulation 

The coil design problem may be solved analytically for some special cases [7J but for other 
geometries the physical problem must be described by a finite number of parameters in 
order to apply numerical methods. The type of parameterisation may chosen to best 
suit the type of coil that is to be designed, ip( r ) can be approximated as a finite weighted 
sum of N basis functions, 

N 

#:)«5>n^n(r), (7) 

n 

and so can the current density by combination of (|7_D with (JBj), 

JV 

J(r)ttJ2i>M*l (8) 



where ^n(r) and j„(r) are the nth stream-function and current density basis functions 
respectively, and ip n are the weights. 

Equation (jSJ) can be incorporated in (j2J) to (jSJ), so that B z (r), W, P and r are 
parameterised as finite summations. 

2.3. Matrix Equations 

This discrete formulation allows matrix equations for each of the physical properties 
to be written. A vector of stream-function weights, t/>, was defined; ip = 
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[ipi, . . . , ip n , . . . , iPn} T (where T represents the transpose operation). Each Cartesian 
component of the current density at a set of points, r s , can be written as matrix equations 

jx = Jxi>i jy = Jylp, jz = Jz^P, (9) 

where j x is a vector that lists values of the x-component of the current density at a set 
of S points, j x = [J x (ri), . . . , J x (r s ), . . . , J x (rs)} and J x in an S x iV matrix. Similar 
matrix equations can be written for the cylindrical coordinate system to give j p , and 

jz- 

A H x N matrix B relates ^ to a vector b of length H containing magnetic field 
values, where b = [Bgfa), . . . , B z (r h ), B z (v H )] T ; 

b = Bip. (10) 

Similarly, each component of the torque vector ()5]) can be written as the inner 
product of if) and a vector, 

T X = T^, Ty = Tylj;, T Z = T^. (11) 

The energy terms ([3]) and (j4|) are quadratic with respect to ■0, 

W = ip T L c ijj (12) 

P = ^ T Rc^ (13) 

where L c and R c are symmetric, N x N matrices of the inductance and resistance of 
the coil surface respectively. In fact, R c = J x J x + JyJ y + J^Jz- 

Three types of parameterisation are used in this work. The first assumes that the 
current-carrying surface is a finite-length cylinder and that ip(r) is a weighted sum of 
truncated sinusoidal functions [HJ[T5]. I n the second approach, surfaces are described 
by flat triangular elements and 0>(r) is a piecewise-linear function [201 123 123 [29] . The 
third approach uses surfaces of revolution about the z-axis and is an axisymmetric BEM 
[2T| [22] . The way in which ^(r) and J(r) are parameterised in each case are given in 
the Appendices. For details of how to calculate the matrices, B, T x , T y , T z , L c , R c 
for each type of parameterisation, the reader is advised to seek the above references. 
Other approaches to discretizing the problem are possible, such as using quadrilateral 
elements, but are not described in the present work. 



2.4- Numerical Problem 

The vector ip is a list of the stream-function weights, ip n , which are the free parameters 
of the coil design problem. The problem including the maximum current density and 
all other terms can be written generally 

mm{UW = /(V) + ae(iP) + (3W^) + + %'WIU}- (14) 

It contains terms to control the residual primary field, f(if>), eddy current field, e(ijj), 
stored magnetic energy, W(i/j), power dissipation, P(ip), and maximum current density, 
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||j(VOIIoo, along with their respective, user-definable weighting factors, a, (3, 7 and 5. 
One, two or three of these parameters are usually set equal to zero to remove them 
from U(ip). For example, 7 = 5 = will result in an actively shielded, torque-balanced 
coil with minimal stored energy. Each term in ( IT4|) possesses a natural scaling from the 
physical constants used in their calculation. Choice of a, (3, 7 and 5 values must balance 
these scalings: for example, a, (3, 7 and 5 are typically in the order 1, 10~ 7 , 10~ 9 and 
10 -10 respectively so that they have a magnitude comparable with the f{ip) term, but 
are dependent on the specific problem. 

The minimization was performed such that if), belonged to the set of stream- 
functions, ^, that exhibit zero net torque (TTTi) ; 

= {ip G R N , T x ifj = & Tyifj = 0}, (15) 

where T x if) and T y if> give the x- and y-components of the torque vector t x and r y , 
respectively. In this work it was assumed that the background magnetic field was 
uniform and oriented parallel to z, and as such r z = 0. It is also possible to balance the 
torque of coils immersed in non-uniform background magnetic fields. 

The / term in ( |T4l) represents the sum-of-squares of the error in the primary 
magnetic field, 

f = ±\\B c iP-b t \\i, 

where || • H2 is the classical £ 2 -norm, B c is a H x iV matrix relating if) to the magnetic 
field values at the H target field points in the ROI ( JTUl) and bt is a vector of length H 
containing the target magnetic field values. 

The term e in f lT^|) represents the sum-of-squares of the magnetic field that the eddy 
currents produce in the ROI. A Heaviside function in coil current was assumed [27] and 
the stream-function of the instantaneously-induced eddy current density, if) e , is linearly 
related to if> by 

ifj e = -L- l M ec if), (16) 

where L e is an N e x N e (N e is the number of basis functions approximating the current 
density on the eddy current surface) self-inductance matrix of the conducting surface 
where eddy currents are induced and M ec is an JV e x iV matrix of the mutual inductance 
between the coil surface and eddy current surface. 

The field produced by the eddy current at the target points was desired to be 
minimal, hence we use the sum-of-squares eddy current field to enforce active magnetic 
shielding. 

e= l -\\B e L~ e 1 M e M\l (17) 

where B e is a H x iV e matrix relating if) e to the eddy current magnetic field values at 
the target points. 
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The stored magnetic energy, W, and power dissipation, P, terms in ( TT41) are 
quadratic with respect to ip and are given by Equations f|T2|) and ffT3l) . respectively 

The maximum current density magnitude in the coil design is written here as the 
£°°-norm of the current density magnitude vector, j. j is a list of length S containing 
the current density magnitude values at each surface point 



2.5. Optimization Algorithm 

Previous methods solved min{?7(^)} by partial differentiation of U(i/j), f^, and 
subsequent matrix inversion of the consequential system of linear equations [28, 29J. 
This cannot be done with (Ti~4"I) since U(ip) contains a non-differentiable £°° term. To 
solve ( !T4l) we used an accelerated descent algorithm of Nesterov [19J on a dual problem 
smoothed using ideas of Moreau-Yosida. Full details of the algorithm will be submitted 
elsewhere. For the purposes of this paper it should be noted that the algorithm requires 
as inputs the smoothing parameter, fi, for the £°°-norm, the number of iterations to 
perform, Q, and an initial guess for the solution, i[> . /i and Q are related by some 
inverse relationship that requires more iterations when less smoothing is applied, but 
will approximate the non-differentiable £°°-norm to a greater degree. Convergence was 
checked by observing the value of the dual cost function as q — > Q. 

The optimization algorithm was coded in Matlab (The Mathworks, Natick, MA) 
and was excecuted on a 64-bit Linux server with Intel (Intel Corporation, Santa Clara, 
CA) Xeon E5430 quad-core CPUs at 2.66GHz. 

2.6. Examples 

The impact of designing coils with the minimax|j| was investigated by three examples. 
These examples are described in the following sections and were chosen to elucidate the 
behavior of the system when designing realistic coils. In the examples outlined below 
the convergence rates and calculation times were recorded. 

Relevant properties of the coil performance were recorded in all cases. The 
efficiency, 77, is the intensity of magnetic field that the coil can generate with 1 Amp and 
is also sometimes referred to as the sensitivity. Inductance, L, resistance, R, minimum 
spacing between wires, w, maximum field error in the ROI, m&x(AB z ) and maximum 
eddy current field in the ROI, ma.x(B ez ), are all recorded. Derived figures-of-merit 
(FOMs) r] 2 /L, rf/R and r/w are independent of the number of contours, N c , used to 
convert ip into wires and are useful for comparing between coils. 




(18) 




(19) 
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2.6.1. Cylindrical X- gradient coils We initially demonstrate minimax \j\ coil design 



with ip(r) parameterised by a sum of sinusoidal basis functions [T4"| [15]. Appendix A 
details the parameterisation of ip{ r ) an d J(r)- The current-carrying surface for these 
examples was assumed to be a finite length cylinder 760 mm in diameter. The region of 
uniformity (ROU) was a 400 mm long, 400 mm diameter cylindrical region positioned 
concentrically inside the current-carrying surface. The target field in the ROI has a 
magnitude that varies linearly in the x-direction; B z (r) oc x. This simple geometry 
was used to investigate some of the more fundamental behaviours of coils designed with 
minimax|j|. In all cases max(AP 2 ) was kept at 5 ± 0.01 % and no torque balancing 
was required because J(r) is forced to be symmetric by the limited parameterisation. 
No active shielding was used, so a = 0. 

The length of the coil, I, was varied between 700 and 2000 mm to observe its 
behavior with respect to standard min(W) and min(P) coils. 

In a second experiment, minimax|j| coils were designed with varying amounts of 
power minimization to investigate their behavior on the continuum from min(P) to 
minimax|j|. This was performed with I = 1400 mm, for A^ = 36 and 200. For A" = 36, 
with reference to (tO) . M' = 4 and JV' = 8. For JV = 200, M' = 10 and JV' = 20. 



2.6.2. Shielded Gradient Coils Short, cylindrical, actively-shielded gradient coils were 
designed with minimax|j|. The dimensions of a coil presented in Reference [30] was 
used in this example. The system was modelled with a triangular BEM [20] which 



approximates ip(r) as piecewise-linear in each triangle as described in Appendix B 



Four different X-gradient coils were designed using different types of minimization, 
all with max(AP 2 ) = 5 ± 0.01% in the ROI and a = 20; min(W / ), min(P), minimax|j| 
and min(P & max|j|) which is some combination of min(P) and minimax|j|. 

A full set of gradient coils comprises X, Y and Z coils, with Y being a 90° rotation of 
the X-gradient. Z-gradient coils (B z (r) cx z) were designed with min(P), minimax|j|and 
min(P and max|j|). It is known that for axi-symmetric geometries and target fields 
(i.e. zonal coils) ip{v) is invariant. All nodes with identical p n and z n were treated 
together and forced to have the same value of ip n . In some way this is equivalent to the 
axisymmetric case described in the next example. 



2.6.3. Shim Coils Designing shim coils not only requires the production of magnetic 
fields that have a different spatial form to gradient coils, but the engineering and 
electronic requirements are also different. The efficiency, rj, is of primary importance 
and higher order shims are considerably less efficient than low order shims. Improving 
rj for the higher order shim coils would be useful to improve correction of geometric 
distortion in MR images induced by B field error and may provide smaller linewidths 
for MR spectroscopy via higher order shimming. Due to the often constrained axial and 
radial space provided for shim coils, wire-spacing can become a problem and limits 77. 
Coils designed with the minimax|j| were studied to see if they could help improve shim 
coil performance. X2-Y2 biplanar shim coils (B z (r) oc x 2 — y 2 ) were designed with 860 
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mm diameter and 500 mm separation [21]. The ROI is a spherical volume of 380 mm 
diameter in which max(AP 2 ) was fixed at 10 ±0.01 %. In both cases, the <fi dependence 



of ip(r) was spectrally decomposed in terms of 7 sinudoids, i.e. M' = 7, see Appendix C 



3. Results 



3.1. General Observations 

Calculation of the system matrices for both the sum-of-sinusoids and the axisymmetric 
BEM took on the order of a few seconds. The time required for the triangular 
BEM system matrix calculations was reduced over previously reported times [29] by 
coding this part in C. It took less than 10 minutes to calculate all the matrices for a 
large problem containing 3200 nodes and 6144 triangles. For a medium-sized problem 
containing 1632 nodes and 3072 triangles the system matrix calculation time was less 
than 2 minutes. 

The optimization algorithm converged in all cases as expected from its deterministic 
nature. The smoothing parameter, /i, for the £°°-norm controls how close the solution 
to the smoothed problem, ijj* is to the true solution, ip* . Practically, we make /i small 
enough so that no observable difference in the solutions is seen for smaller /i. The time 
required to find a solution close to ip* varies widely and is dependent on the problem. 
The number of iterations, Q, that the algorithm required to converge is inversely related 
to \i. For a typical problem tackled in this work, 1 x 10~ 16 </i<lx 10~ 14 resulted in 
indistinguishable solutions which typically required 10, 000 < Q < 200, 000. 



3.2. Cylindrical X-gradient coils 

The time required to find the solution to a small problem with /i = 1 x 10~ 15 , Q = 20, 000 
and the number of free variables, N = 36 was approximately 23 seconds. Figured] shows 
how the FOMs for a) stored energy, b) power dissipation and c) wire spacing varied for 
min(W), min(P) and minimax|j| X-gradient coils as the length of the coil surface varied 
with max(AP z ) = 5 ± 0.01%. Figure [2] a) shows the values of 5 that were used with 
varying 7 to maintain max(AP 2 ) = 5 ± 0.01% in the ROU. The variation of the two 
relevant FOMs with 7 are shown in Figure [2] b) and c). Figure [3] shows one quadrant 
of the wire paths for four coils designed with a) min(iy), b) min(P), c) minimax|j| 
with N = 200 and d) minimax|j| with N=36. Wire positions are unwrapped from 
their cylindrical shape onto a flat z-a<fi plane. As with all coils presented in this paper, 
connections must be made during construction from each loop to its neighbour to ensure 
current flow throughout the coil. The location of these coils are marked for reference on 
Figure [2] b) and c) where ®, © and <D to the coils in Figure [3] b), c) and d) respectively. 



3.3. Shielded Gradient Coils 

One quadrant of the wire paths for minfW), min(P), minimax|j|, and min(P & max|j|) 
shielded X-gradient coils are shown in Figure HI The left-hand side show the primary 




Figure 1. Variation of a) ij 2 /L, b) r/ 2 /R and c) rjw with length-to-diameter ratio, l/d 
for unshielded X-gradient coils designed with a sum-of-sinusoids parametrised stream 
function. 



coils and the right-hand side are their active magnetic shields. Their performance 
characteristics are given in Table [U Due to the high number of free-parameters, 
N = 1985, it took approximately 480 minutes to perform Q = 210,000 iterations to 
obtain a well converged solution. 

The stream-functions of the current densities of the three Z-gradient coils with 
min(P), minimax|j|, and min(P & max|j|) are shown in Fig. [5] The FOMs are given 
in Table [2j The calculation time required to find the optimal solution was dramatically 
reduced for the Z-gradient coils by reducing the number of free variables. The time 
required for Q = 80, 000 iterations was 10 minutes for a well converged solution. 
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Figure 2. a) 5 values required in combination with varying 7 in order to maintain 
field error of 5 ± 0.01 % with /? = 0. b) the resulting r/ 2 /R and c) r/w of the coils for 
32 and 200 sinusoidal basis functions. Data points labelled ®, ® and © correspond to 
the coils shown in Figure [3] b), c) and d) respectively. 



3-4- Shim Coils 

The wire-paths of one plane of a minimax|j'| X2-Y2 biplanar shim coil are shown in 
Fig. [6] b) next to an equivalent min(P) coil. Given a 4 mm wire spacing limit for 
construction, the maximum achievable 77 were 73.8 and 102.0 mTm" 2 A _1 for the min(P) 
and minimax|j| coils using iV c = 17 and 20, respectively. 
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-0.4 -0.2 0.2 0.4 -0.4 -0.2 0.2 0.4 



C) d) 




aip (m) acp (m) 

Figure 3. One quadrant each of the wire paths for the sum-of sinusoids, unshielded, 
X-gradient coils in the cases of a) min(T / l /r ), b) min(P), c) minimax|j| with 200 
sinusoids and d) minimax |j| with 32 sinusoids. Red wires indicate reversed current 
flow with respect to blue and only 12 contours of the stream function are shown for 
clarity. 

4. Discussion 

This paper reports a method to directly minimize the maximum current density for the 
magnetic field synthesis or coil design problem. It focuses on the design of gradient 
and shim coils for MRI applications, but can be considered as a general approach to 
coil design; it may prove useful outside the realm of gradient and shim coil design. 
Superconducting magnet design (e.g. see Ref. [32]) is one such application that merits 
some comment. Although no experiments have yet been performed on magnet design 
using this minimax optimization, it might be employed to design magnets with reduced 
peak current density and/or peak magnetic field in the conductors, for example. For 
practical designs, £ -norm minimisation of the current density may be incorporated to 
yield low peak, yet sparse current density designs. 

Previous gradient coil design methods have very effectively minimized the stored 
energy [9] and power dissipation [16] subject to the production of magnetic fields of a 
prescribed accuracy. Other approaches that lower the maximum current density have 
been presented [211 231 [31] , but none can be shown to be optimal in terms of minimax 
current density. We used the adaptive regularisation technique [23] to design shielded X- 
gradient coils shown in Fig. H] and found that it could achieve a minimum wire spacing 
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Figure 4. Wire paths for the triangular BEM, actively-shielded, X-gradient coils in 
the cases of a) and b) min(W / ), c) and d) min(P), e) and f) minimax|j| and g) and 
h) min(P & max|j|). The active magnetic screens appear on the right for the primary 
coils on the left. Red wires indicate reversed current flow with respect to blue and only 
12 contours of the stream function are shown for clarity. 



of 8.7 mm. The mmimax|j| algorithm achieved 9.3 mm wire spacing indicating that 
adaptive regularisation works well, but cannot maximally spread the wires. The reason 
for the difficulty in achieving truly minimax current density coils is that such a term is 
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Table 1. Properties of the triangular BEM designed, actively-shielded X-gradient 
coils. Input parameters a, /3, 7, S and coil properties including efficiency, 77, maximum 
field error in the ROU, max(Ai? z ), maximum eddy current field in the ROU, max(_B e2 ), 
inductance, L, resistance, i?, minimum wire spacing, w and figure-of-merit values, 
rj 2 /L, rj 2 /R, r]W. 
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Figure 5. Stream-functions for the actively-shielded, Z-gradient coils designed with 
min(P), minimax I j' I and min(P & max|j|). The shield stream functions are greater 
than zero and the primary stream function are less than zero. 

non-differentiable with respect to the solution variables. It is surely possible to insert 
such a non-differentiable term into a stochastic optimization technique such as a genetic 
algorithm [33] or simulated annealing [34], but it is expected that such methods would 
require very long computing times and converge to a solution that is not necessary the 
global one. Here, the maximum operation is expressed as the infinity norm (also known 
as the uniform or Chebychev norm), || • H^, smoothed and converted to its norm-dual, 
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Table 2. Properties of the triangular BEM designed, actively-shielded Z-gradicnt 
coils. Input parameters a, /3, 7, 5 and figure-of-merit values, r] 2 /L, r] 2 /R, rjw are 
given. 
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Figure 6. One plane of the wire paths for the axisymmetric BEM designed biplanar 
X2-Y2 shim coils in the cases of a) min(P) and b) minimax \j\. Red wires indicate 
reversed current flow with respect to blue. 

the £ 1 -norm. 

The time required for this algorithm to converge varies widely on the size of the 
(S x N) current density matrices, J x , J y and J z . For the Turner- Carlson approach 
[H El HH EE] with iV = 32 and S = 441 convergence was obtained in 23 seconds. 
However, for iV = 1985 and S = 4096 with the triangular BEM it took 480 minutes. It 
should be noted that for 6 = 0, the solution is obtained in less than a second as just 
one matrix inversion is needed (281 129] . This illustrates the need to take into account 
any symmetry that might be present in the system to lower the number of free variables 
and speed up the calculation. For a defined maximum field error, the design process 
needs to be repeated in order to obtain the ideal trade-off parameter, 5. It is hoped that 
the amount of user input and computational burden can be reduced by describing the 
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problem as a constrained one in which maximum field error is a user-definable constraint. 

The length of an unshielded X gradient coil was varied such that the length-to- 
diameter ratio, l/d, ranged from 0.92 to 2.63. The performance of the coils designed with 
min(VF), min(P) and minimax|j| were evaluated and Figure [1] shows the dependencies 
of rj 2 /L, r] 2 /R and rjw on coil length. rf /L characterises the power requirements of the 
driving amplifier and rf /R characterises the total amount of heat generated by the coil 
where higher values indicate better performance in both cases, rjw on the other hand 
characterises the maximum field strength that can be obtained irrespective of inductance 
or resistance for a given minimum wire spacing. Several interesting behaviours are 
apparent from studying the data in Figure [TJ First, when designing a coil with min(VF) 
it will have the highest rf /L value. Likewise, a min(P) coil will have the highest rj 2 / R 
value and minimax|j| coils will have the highest rjw. This is to be expected. min(W) 
and min(P) designed coils have very similar rf / L values and slightly different rf / R 
values, with the minimax|j| coils possessing lower values of these two FOMs. In fact, 
for very long coils (l/d > 2.2), the value of rf/L and rf /R actually decreases for the 
minimax|j| coils. However, the minimax|j| coils possess a rjw value that is considerably 
larger than that of the min(VT) and min(P) designed coils. Figure Q] c) shows that for 
minfW) and min(P) coils with l/d > 2 the value of r]w becomes flat. This happens 
when the region of max|j| occurs approximately at the end of the ROI and not at the 
end of the coil. This indicates that the length of the coil surface is no longer restricting 
the maximum achievable field strength when l/d > 2. rjw appears to be tending to a 
particular value for long minimax|j| coils that is approximately 1.5 to 2 times larger 
than the other coils. Unlike a previous approach [23], minimax current density coils 
are dramatically different from their Tikhonov regularized counterparts even for long 
cylindrical coils. 

It is known that a unique solution is found when ill-posed problems are solved with 
Tikhonov regularisation, which is the case for min(W / ) and min(P). It is not known if 
a unique solution results from including the minimax|j| term in the functional. It is 
suspected by the authors that there is no unique solution, but more theoretical analysis 
is required to establish this. Figure |2] shows the behavior when /3 = and both 7 and 
5 are finite, i.e. as min(P) is traded for minimax|j| in the optimisation. The 5 value 
required to maintain a constant field error is inversely related to 7, as expected. 7 and 
5 values are similar for N = 32 and 200 sinusoidal basis functions. Figure |2] b) shows 
the variation of the power FOM, rf /R, as this trade-off happens. It is evident from 
these data that by adding a small amount of 7 to the minimax|j'| coil, a sharp increase 
in rf /R can be effected at the expense of very little decrease in r]w. When N = 32 more 
smoothness is enforced by the basis functions and P is limited. Hence rf /R is lower 
when more basis functions are used, indicated by © and ® on Figure [2] b), but rjw is 
also limited. This difference is also evident in Figures |3]c) and d). Conversely, Figure 
[2] c) shows that by adding a small amount of 5 to the min(P) coil, a large increase in 
7]w can be achieved with only a small change in the power dissipation of the coil. It is 
not surprising to observe that both min(P) coils with N = 32 and 200 are essentially 
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the same since min(P) coils favour low spatial frequencies in ?/>(r). It is clearly possible 
to choose any solution on the continuum from min(P) to minimax|j|. Although not 
presented in this paper, it is also possible to trade min(W) with minimax|j| along a 
similar continuum. 

From Figure |3]it can be seen that the min(W) coil possesses an area with the highest 
current density at the ends of the primary with the wires of the power minimized coil 
being more spread, as expected [16]. ip(r) conforms to the usual cos (ft dependence for 
the min(W / ) and min(P) coils despite no such constraint. However, this leads to regions 
of higher current density at <fi = which is optimally dispersed in the minimax|j| 
coils. Deviation from the cos</> behavior at the ends of the coils means that spherical 
harmonic fields of higher degree will be introduced in the ROI. These high degrees are 
then cancelled by ip(r) variations closer to the ROI. 

Incorporating active magnetic shielding (26] into the functional is a simple matter 
since it can be written in a similar form as the target field term. Actively-shielded 
X-gradient coils were designed using the same geometry as appears in Reference [30J. 
Figure H] shows one quadrant each of the primary and shield wire paths of the min(W), 
min(P), minimax|j| and min(P & max|j|). It is interesting to note that the shield coil 
for the purely minimax|j| coil has unnecessary current density with many reversed turns, 
Figure H]f). This results from the fact that there is no penalty for extra current density 
when /3 and 7 are zero. By incorporating a small amount of 7, this impractical design 
is very effectively converted into a highly practical design with smooth wire paths, low 
resistance and very well spread wires, Figure H]h). It would also be simple to have 
different 7 and S values for primary and shield coils. The wire paths in Figure H]e) show 
the tendency of minimax|j| in an extreme case where right-angular corners appear in 
the design. 

Similar, but less pronounced effects were observed from the results of the Z-gradient 
coil of identical geometry. Figure |5] shows the stream functions along the z-direction for 
the coils designed with max(AP 2 ) = 5 ± 0.01 %. The magnitude of the current density 
(in this case the steepness of the slope of the stream-function) is almost the same in 
all parts of the coil. This again leads to unnecessarily large amounts of current density 
on the shield coil, which is easily removed by the addition of a small value of 7. The 
combined min(P & max|j|) Z-gradient coil exhibits large wire spacing and marginally 
increased resistance when compared to the min(P) coil. The problems associated with 
high current densities are less severe when compared to those of X-gradient coils, but 
this approach may be more useful for zonal shim coils of higher order. 

In a final example, an axisymmetric BEM was used to design biplanar X2-Y2 shim 
coils. Rotational symmetry of the system about the z-axis is assumed. ip{ r ) is spectrally 
decomposed in and spatially in p and z. Qualitatively, it can be seen from Figure 
[6] that the minimax|j| coil used all the space provided, whereas the min(P) coil forced 
wires to be very smooth. For a fixed max(AP 2 ) = 10 ± 0.01 % and w > 4 mm, n is 38% 
higher for the minimax|j| coil. The construction of such coils may be made slightly more 
complex by the additional loops in the design. Again, a combined min(P & max|j|) coil 
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might provide a good balance between simplicity and efficiency. 

5. Conclusion 

It has been shown that the magnetostatic field synthesis problem can be solved for coils 
with minimax current density. The problem was solved in the present study with three 
different boundary element methods to illustrate the generality of the approach. It 
can therefore be used to synthesize any physically realistic magnetic field with currents 
flowing on arbitrary surfaces. Alternatively, the time needed to solve the problem can 
be dramatically reduced by assuming some degree of symmetry. Coils with minimax 
current density possessed increased resistance and inductance for the same field error 
and in some cases had high current densities in regions known to naturally require only 
low current density. More practical coils were obtained with a mixture of power and 
maximum current density minimization. Such coils are characterized by low inductance 
and resistance, but also a large spacing between all the wires of the coil. This spreading 
of wires may be used to increase the efficiency of the coil by permitting extra turns to 
be added, reduce turn-to-turn eddy current effects, reduce localized heating in the coil, 
and design coils that are easier to manufacture and/or are extremely short. Moreover, 
it allows a coil designer to explore a new range of optimal solutions to the field synthesis 
problem. The resultant coils must be judged by figures-of- merit appropriate to the 
desired characteristics of the coil, since field accuracy, gradient efficiency, stored energy, 
power dissipation and maximum current density may all be traded for each other. 
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Appendix A. Sum-of-Sinusoids 

In this appendix we present for completeness the formulation for calculating the matrices 
in Equations (JSJ) to ffl3|) necessary for implementing the minimax current density 
algorithm with sinusoidal stream-function basis functions. In this case, the coil surface 
is assumed to be a finite-length, /, cylinder of radius a with its axis of symmetry 
oriented in the z direction. The stream-function of the current density, is spectrally 
decomposed and assumed to be a finite weighted sum of truncated sinusoidal functions 
in z and (ft; 



M' 



N' 




(A.l) 
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where 



/27m' \ 

( — J — z j cos ((2m' — 1)0) if \z\ < ^ 



sin 

vw«'(&*) = <( v ^ y 1 1 ' . (a.2) 

if \z\ > | 

The prime indicates the difference between n used in the main algorithm and the 
order, n', and degree, m', of the sinusoid. Equation ([7]) is obtained by a reordering 
of the weights \ m 'n' as ip n . It restricts the magnetic field to be antisymmetric in the 
x-direction and symmetric in the z-direction. Extra basis functions that are symmetric 
in x and antisymmetric in z can be included in order to remove the inherent field 
symmetry enforcement pj], [12j [13] , but this is not required in the present study since 
we are designing an X-gradient coil and know that ^(r) = at \<p\ = ~. Due to this 
enforced symmetry it is known that the net torque experienced by the coil is zero. 

The current density on the coil surface has J^- and J 2 -components that are 

J^cp, z) = = A m / n /^p- cos [—J- Z J cos (( 2m ' - 1)0) ( A - 3 ) 

ml n' 
m! n' 

and are equivalent to (jHJ) and fl2]). 

Appendix B. Triangular Boundary Elements 

A surface can be meshed as a series of / triangular elements with N nodes at the 
corners of the triangles [20]. In this case, ip{v) is piecewise-linear in each triangle and 
the stream- function values at the node positions, if) n , define the whole stream-function; 

N I 

if>(r) = J2^I2^i(r) (B.l) 



n=l i=l 



where 



Mr) = 1 - (r " ra) ' dm (B.2) 

I £*ni | 

if r is a point in triangle i and n is a node of that triangle. ip n i(r) = otherwise. It 
is possible to use higher order shape functions over the triangle so long as they form a 
divergence-free basis [35] . 

The current density on the surface is found from (IB.lj) . (IB. 21) and (jSJ) yielding ([8]) 

and 

J>) = £v ra (r) = ]T|^ (B-3) 
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if r is a point in triangle % and n is a node of that triangle. v ni (r) = otherwise. Ai is 
the area of the triangle i and e ni is the vector that describes the edge of the zth triangle 
opposite the nth node. This demonstrates that the current density is uniform over each 
element of the mesh and so that there needs to be one current density sample for each 
/ triangles in order to fully characterise it. Therefore s becomes i and S is equal to the 
number of triangles, I. 



Appendix C. Axi-symmetric Boundary Elements 



The axisymmetric BEM is used for coil supports that can be described by surfaces of 
revolution about the z-axis. Each "node", n', of this surface is in fact a circle in the 
xy-p\&ne and defined by its radius, p n > and axial position z n >. There may be a conical 
element, i, either side of each node, labelled + and — . A local coordinate, ((p,z), is 
defined for each element that is at one end and 1 at the other. Positions on these two 
conical surfaces are 
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The surface and therefore ip(r) are parameterised in p and z. ip{r) is decomposed 
spectrally parameterised in the ^-direction as a sum of sinusoids. In section §2.6.31 
an X2-Y2 shim coil is designed that has a target magnetic field with 2-fold rotational 
symmetry about z. Therefore, ip(r) is restricted in the 0-direction to take the form 
cos (2(2m' — 1)0). As described in Appendix A the basis function weights, A m v can 
be reordered to comply with the vector arrangement, ip, in the main algorithm. 

M' N' I 

VKC, 0) = Yl Xm ' n ' 2 ^m'n'i(C) cos (2(2m' - 1)0) (C.2) 

m' n' i 

where 

VWn'i(C) = (1 - (C3) 

if i is on the positive side of n', for < ( < 1, 

Vwi(C) = C (C4) 

if i is on the negative side of n', for < C < 1 and ip m 'n'i(C) = otherwise. 
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The stream-function §6§ is applied to obtain, after considerable amounts of algebra, 
the discretised current density 



where 



M AT 



J ( r ) = ^ m ' n ' v ^'^( r ) 



(C.5) 



V m 'n'i(r) 



+ 



cos(2(2m' — 1)0) sin0 

(Pn+1 ~ Pnf + (Z n +1 - Z n ) 2 

m'jl - C)(Pn+i - Pn) sin(2(2m' - 1)0) cos0 

^ (p n+ i - p n ) 2 + 0„+i - Z„) 2 (p n + C(Pn+l - Pn)) 

— cos(2(2m' — 1)0) cos0 

\/ (Pn+l - Pn) 2 + (Z n+ i - Z n ) 2 

| m'(l - C)(p»+i - Pn) sin(2(2m / - 1)0) sin0 

-sj {p n+ i - p n ) 2 + (Z n+ i - Z n ) 2 (p n + C(Pn+l - Pn)) 

[ ~ (z n+1 - z n ) sin(2(2m' - 1)0) 

\J (Pn+1 - Pn) 2 + {Z n +1 ~ Z n ) 2 (p n + C(Pn+l ~ Pn)) 

if i is on the positive side of n', for < ( < 1 and < < 2n, 



z (C.6) 
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if i is on the negative side of n', for < ( < 1 and < < 2n and v m /„/j(r) = 
otherwise. 
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